function [min_male, max_male, min_female, max_female, problem] = get_povrate(couple,singlemale,singlefemale,leisurem,leisuref,hworkm,hworkf,wagem,wagef,q,Q,nonlabor,mutually_acceptable,divorce_costsM,divorce_costsF,divorce_costsMF,hcat,pline,text_solver)

% ================
% define constants
% ================

T = 112;                                              % total potential working hours
H = length(leisurem);                                 % total number of households in the data
alpha = 0.4;                                          % bound for non-labor income

% =========================
% define decision variables
% =========================

yalmip('clear')

qM = sdpvar(H,1);                                     % private consumption of Hicksian good male
qF = sdpvar(H,1);                                     % private consumption of Hicksian good female

pM = sdpvar(H,H,'full');                              % personalized public good price male
pF = sdpvar(H,H,'full');                              % personalized public good price female

pmhworkM = sdpvar(H,H,'full');                        % personalized male's hwork price male
pmhworkF = sdpvar(H,H,'full');                        % personalized male's hwork price female

pfhworkM = sdpvar(H,H,'full');                        % personalized female's hwork price male
pfhworkF = sdpvar(H,H,'full');                        % personalized female's hwork price female

wageM = sdpvar(H,1);                                  % shadow wage male
wageF = sdpvar(H,1);                                  % shadow wage female

nl_incomeM = sdpvar(H,1);                             % non-labor income male
nl_incomeF = sdpvar(H,1);                             % non-labor income female

sharem = sdpvar(H,1);                                 % consumption based sharing-rule male
sharef = sdpvar(H,1);                                 % consumption based sharing-rule female

povertym = binvar(H,1);
povertyf = binvar(H,1);


%%
% ==================
% define constraints
% ==================

%%
% ==================
% define constraints
% ==================

Constraints = [];

% restriction on individual private consumption
for i=1:H
    Constraints = [Constraints, q(i) == qM(i) + qF(i), qM(i) >=0, qF(i) >= 0];

    if singlefemale(i) == 1
        Constraints = [Constraints,qM(i) == 0];
    elseif singlemale(i) == 1
        Constraints = [Constraints,qF(i) == 0];
    end

end


% shadow price should be non-negative and equal to observed wage if the individual is employed
for i=1:H

    Constraints = [Constraints,wageF(i) >= 0, wageM(i) >= 0];

    if ~isnan(wagem(i))
        Constraints = [Constraints, wageM(i) == wagem(i)];
    end
    if ~isnan(wagef(i))
        Constraints = [Constraints, wageF(i) == wagef(i)];
    end

end


% restriction on public prices
for i=1:H
    for j=1:H

        Constraints = [Constraints, 1 == pM(i,j) + pF(i,j), pM(i,j) >=0, pF(i,j) >= 0];
        Constraints = [Constraints, wageM(i) == pmhworkM(i,j) + pmhworkF(i,j), pmhworkM(i,j) >=0, pmhworkF(i,j) >= 0];
        Constraints = [Constraints, wageF(j) == pfhworkM(i,j) + pfhworkF(i,j), pfhworkM(i,j) >=0, pfhworkF(i,j) >= 0];

        if (couple(i) == 1 || singlemale(i) == 1) && singlemale(j) == 1
            Constraints = [Constraints, 1 == pM(i,j), wageM(i) == pmhworkM(i,j), wageF(j) == pfhworkM(i,j)];
        end

        if (couple(j) == 1 || singlefemale(j) == 1) && singlefemale(i) == 1
            Constraints = [Constraints, 1 == pF(i,j), wageM(i) == pmhworkF(i,j), wageF(j) == pfhworkF(i,j) ];
        end
    end

end


% feasibility restrictions on the individual nonlabor income (for each matched pair)
% individual nonlabor incomes after divorce must lie between 40% and 60% of total nonlabor income under marriage
for i=1:H

    Constraints = [Constraints, nonlabor(i) == nl_incomeM(i) + nl_incomeF(i)];

    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints, alpha*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= (1-alpha)*nonlabor(i), alpha*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= (1-alpha)*nonlabor(i)];
        else
            Constraints = [Constraints,(1-alpha)*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= alpha*nonlabor(i), (1-alpha)*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= alpha*nonlabor(i)];
        end
    elseif singlefemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
    end

end


% individual rationality
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, wageM(i)*T + nl_incomeM(i) - divorce_costsM(i) <= qM(i) + wageM(i)*leisurem(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
        Constraints = [Constraints, wageF(i)*T + nl_incomeF(i) - divorce_costsF(i) <= qF(i) + wageF(i)*leisuref(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
    end
end


% no blocking pairs
for i = 1:H
    for j = 1:H
        if mutually_acceptable(i,j) == 1
            Constraints = [Constraints, (wageM(i) + wageF(j))*T + nl_incomeM(i) + nl_incomeF(j) - divorce_costsMF(i,j) <= ...
                qM(i) + qF(j) + wageM(i)*leisurem(i) + wageF(j)*leisuref(j) + pM(i,j)*Q(i) + pF(i,j)*Q(j) + ...
                pmhworkM(i,j)*hworkm(i) + pmhworkF(i,j)*hworkm(j) + pfhworkM(i,j)*hworkf(i) + pfhworkF(i,j)*hworkf(j)];
        end

    end
end


% define sharing rule
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, sharem(i) == qM(i) + Q(i),...
            sharef(i) == qF(i) + Q(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, sharem(i) == q(i) + Q(i)];
    elseif singlefemale(i) == 1
        Constraints = [Constraints, sharef(i) == qM(i) + Q(i)];
    end
end


% define poverty indicator
for i = 1:H
    Constraints = [Constraints, 1 - sharem(i)/pline + 0.001 <= 10000*povertym(i),...
        10000*(povertym(i) - 1) <=  1 - sharem(i)/pline,...
        1 - sharef(i)/pline + 0.001 <= 10000*povertyf(i),...
        10000*(povertyf(i) - 1) <=  1 - sharef(i)/pline];
end


%%
% =================
% Solve the problem
% =================

min_male = zeros(16,1);
max_male = zeros(16,1);
min_female = zeros(16,1);
max_female = zeros(16,1);

problem = 0;

options = sdpsettings('verbose',0,'solver',text_solver);

% males
j = 1;
for empym = 1:2
    for edum = 1:2
        for empyf = 1:2
            for eduf = 1:2

                k = empym*1000 + edum*100 +  empyf*10 + eduf;
                index = (couple == 1 & hcat == k);

                if sum(index) > 0 && problem == 0

                    obj1 = sum(povertym.*index)/sum(index);
                    obj2 = sum(povertyf.*index)/sum(index);


                    % minimize males' share
                    sol = optimize(Constraints,obj1,options);
                    if sol.problem ~= 0
                        sol.info;
                        yalmiperror(sol.problem);
                        problem = 1;
                    else
                        min_male(j) = value(obj1);
                    end

                    % maximize males' share
                    sol = optimize(Constraints,-obj1,options);
                    if sol.problem ~= 0
                        sol.info;
                        yalmiperror(sol.problem);
                        problem = 1;
                    else
                        max_male(j) = value(obj1);
                    end

                    % minimize females' share
                    sol = optimize(Constraints,obj2,options);
                    if sol.problem ~= 0
                        sol.info;
                        yalmiperror(sol.problem);
                        problem = 1;
                    else
                        min_female(j) = value(obj2);
                    end

                    % maximize females' share
                    sol = optimize(Constraints,-obj2,options);
                    if sol.problem ~= 0
                        sol.info;
                        yalmiperror(sol.problem);
                        problem = 1;
                    else
                        max_female(j) = value(obj2);
                    end

                else
                    min_male(j) = NaN;
                    max_male(j) = NaN;
                    min_female(j) = NaN;
                    max_female(j) = NaN;
                end
                j = j + 1;
            end
        end
    end
end



return


